function [elements,processInfo] = CalculateMaterialResponse...
    (elements,model,processInfo)
ne = model.countElements();
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    ue = f.getUHatLocal();
    ng = e.getNumGaussPoint();
    B = zeros(3,8);
    C = elements(i).constitutiveMatrix;
    if elements(i).isCrack == 0
        for j=1:ng
            xi = e.getGaussPoint(j-1);
            B = e.computeB(f,B,xi);
            %compute stress
            elements(i).stress(:,j) = C*B*ue;
        end
    else
        %% update Xi
        due = ue - elements(i).ue_old;
        dxi = elements(i).Kxixi\(elements(i).renh-...
            elements(i).Kxid*due);
        elements(i).Xi = elements(i).Xi+dxi;
        
        %% update stress at each Gauss point
        for j=1:ng
            xi = e.getGaussPoint(j-1);
            B = e.computeB(f,B,xi);
            %compute strain
            Gc = computeGc(elements(i),B);
            epsilon = B*ue + Gc*elements(i).Xi;
            %compute stress
            elements(i).stress(:,j) = C*epsilon;
        end
        
        %% compute and update renh
        elements(i) = computeRenh(e,elements(i));
        
    end
    %% save ue
    elements(i).ue_old = ue;
end
end
